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Preface 

In  March,  1078,  the  Air  Force  Weapons  Laboratory  proposed  a thesis 
topic  which  involved  calculating  the  electric  field  inside  a spacecraft 
which  had  been  exposed  to  a field  of  incoming  current  density.  This 
problem  interested  me  for  two  reasons.  First,  it  has  a weapons  appli- 
cation, and  the  results  of  this  project  would  be  useful  to  the  sponsor- 
ing office  at  AFWL.  Secondly,  working  on  this  type  of  problem  would 
strengthen  my  background  in  electromagnetism. 

I would  like  to  express  my  sincere  gratitude  to  my  advisor, 

Dr.  Donn  G.  Shankland.  Without  his  guidance  this  project  would  not 
have  been  possible.  I would  also  like  to  thank  Dr.  David  A.  Lee  for  his 
assistance  with  variational  calculus,  and  for  his  explanation  of  logarithmic 
potentials.  Special  thanks  are  also  due  to  1LT  Dave  Hardin,  who  helped 
clarify  difficulties  encountered  with  natural  boundary  conditions,  and  to 
Dr.  John  Jones,  who  reviewed  and  critiqued  my  analytical  solution. 

Richard  A.  Passow 


11 


1 


s 

f 


i 


Contents 


List  of  Figures v 

List  of  Symbols vi 

Abstract vii 

I.  Introduction  1 

II.  Peve lopmont  of  a Theoretical  Model  3 

Physical  Analysis  3 

Method  of  Solution  4 

Boundary  Conditions  6 

III.  Solution  by  Separation  of  Variables  10 

Solution  in  Region  I 10 

Solution  in  Region  II 10 

£-Field  in  Region  II 11 

Test  Solution  in  Region  I 11 

IV.  Variational  Calculus  Approach  13 

Solution  Procedure  13 

Matrix  Form  of  the  Double  Integral 15 

Matrix  Form  of  the  Single  Integral 19 

Summary 21 

V.  Solution  Using  Green's  Functions  23 

Reduction  to  Integral  Equation  23 

Solution  in  Region  II 26 

Calculation  of  the  l?-Field 27 

Summary 27 

VI.  Representative  Values  of  the  Magnitude  of  the  E-Field  ...  28 

VII.  Conclusions  and  Recommendations  31 


Bibliography 
Appendix  A: 
Appendix  B: 


Solution  by  Separation  ot  Variables  

Fourier  Expansion  of  the  Boundary  Conditions  . . . . 

iii 


33 

34 


41 


Contents 


Page 


Appendix  C:  Computer  Program  to  Calculate  the  Magnitude  of  the 

li-Field  ,J 45 

Appendix  D:  Derivation  of  the  Minimization  Functional  49 

Appendix  E:  Matrix  Form  of  the  Variational  Equation  53 

Appendix  F:  Importance  of  Natural  Boundary  Conditions  when 

Minimizing  a Functional  58 

Appendix  G:  Variational  Computer  Program  63 

Appendix  H:  Derivation  of  the  Green's  Function  and  Iteration 

Formulas 77 

Appendix  I:  Derivation  of  Integral  Equation  for  <f(x) 86 

Appendix  J:  Green's  Function  Computer  Program  90 

Vita 94 


iv 

J 


Figure 

2-1  Schematic  Drawing  of  the  Theoretical  Model 


Page 


5 

4-1  Scalloping  Phenomena  When  Continuity  of  the  1st  Derivative 


Is  Not  Enforced  Across  All  Mesh  Points 18 

4-2  Plot  of  <£>/cos  $ vs  Theta  Along  the  Outer  Boundary.  Con- 
tinuity of  the  First  Derivative  in  the  Angular  Direction 
Is  Not  Enforced 20 

4-3  Comparison  of  the  Analytic  Solution  with  the  Variational 

Solution  Along  the  Line  0 = 0 22 

6-1  Plot  of  | E j max  vs.  Cylinder  Thickness  for  an  Aluminum 

Cylinder  with  a Radius  of  ,5m 29 

6-2  Plot  of  [E|  vs.  Angular  Displacement  for  an  Aluminum 

Cylinder  with  a Radius  of  . 5m 30 

E-l  Sample  Mesh  for  Calculating  ij>  in  the  Cylindrical  Shell  . . 54 

F-l  Potentials  along  the  Line  9 = 0 as  Solved  by  Imposing 
the  Boundary  Conditions  with  Lagrange  Multipliers  and 
Ignoring  the  Natural  Boundary  Conditions  61 

F-2  Solution  of  a 1-D  Problem  Using  Lagrange  Multipliers 

and  Natural  Boundary  Conditions  62 

H-l  Schematic  for  Calculating  the  Singular  Portion  of 

Eqn  (11-17) 81 

11-2  Schematic  Diagram  for  Calculating  <j>(x)  on  the  Inner 

Boundary  83 

1-1  Domain  for  the  Derivation  of  the  Integral  Equation  ....  86 


I 


v 


A 


A 


B 

5 

— V 

E 


Si-S2 


I 


J 

L 

q 

•p 

s 


o 


1 


,o 


2 


* 


List  of  Symbols 


Represents  the  area  bounded  by  "S"  when  used  as  a limit  of 
integration. 

Double  underline  denotes  a matrix. 

Single  underline  denotes  a column  vector. 

Represents  a delta  function  or  a variance. 

Electric  field. 

Arbitrary  constants. 

Incoming  current  density. 

Current  density. 

Brackets  represent  a functional. 

Lower  triangular  matrix. 

represents  the  transpose  of  a matrix. 

Potential . 

Represents  a functional  or  an  inhomogeneous  source  term.  When 
used  as  a limit  of  integration  it  represents  the  boundary  of  a 
two-dimensional  surface. 

Conductivity. 

Arbitrary  constants. 


vi 


r 1 

Abstract 


If  a spacecraft  is  exposed  to  a steady  stream  of  current  density,  the 
charge  of  which  is  deposited  on  the  surface  of  a cylindrical,  conducting 
spacecraft,  internal  electromagnetic  fields  arc  generated.  If  the  internal 
fields  are  of  sufficient  strength,  undesirable  electronic  noise  or  damage 
may  result.  This  thesis  presents  three  approaches  for  calculating  the  induced 
E-Field:  separation  of  variables,  variational  calculus,  and  the  use  of  Green's 

funct ions . 

The  spacecraft  is  modeled  as  a hollow,  infinite  cylinder.  The  fields 
are  calculated  for  the  case  in  which  the  incoming  current  is  incident  per- 
pendicularly to  the  longitudinal  axis  of  the  cylinder.  Basic  electro-static 
theory  reveals  that  the  governing  equation  for  the  potential  is  Laplace's 
Equation,  subject  to  Neumann  boundary  conditions.  This  equation  is  first 
solved  by  separation  of  variables.  The  E-Field  predicted  for  representive 
values  of  incoming  current  density  are  on  the  order  of  10  7 volts/meter. 

Several  pitfalls  encountered  with  the  variational  approach  are  explained. 
These  include  the  importance  of  natural  boundary  conditions,  ensuring  continu- 
ity of  the  first  derivative  across  all  mesh  points,  and  difficulties  encountered 
in  trying  to  reduce  the  size  of  the  matrix  through  "decoupling".  The  results 
of  this  section  show  that  the  variational  method  tends  to  approximate  the 
analytic  solution  as  the  mesh  becomes  finer. 

The  Green's  function  approximations  of  the  potential  distribution  were 
consistent  with  analytic  results  to  at  least  two  significant  figures. 

vii 
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1 Introduction 


This  thesis  concerns  the  calculation  of  the  electric  field  induced 
inside  a spacecraft  (modeled  as  a cylinder),  which  has  been  irradiated 
by  a steady  stream  of  current  density.  The  source  of  the  incoming 
current  is  not  a factor  under  consideration.  This  problem  was  proposed 
by  the  Air  Force  Weapons  Laboratory,  and  this  paper  provides  an  examina- 
tion of  several  approaches  to  the  problem's  solution. 

This  project  is  limited  to  the  study  of  a non-rotating  cylinder 
which  is  exposed  to  a steady  stream  of  current,  incident  perpendicularly 
to  the  cylinder's  longitudinal  axis.  Also,  the  cylinder  is  considered 
to  be  infinite.  This  reduces  the  problem  to  two  dimensions.  Thus,  the 
answers  calculated  will  not  give  the  exact  fields  inside  an  actual 
spacecraft.  However,  the  answers  are  useful  as  "first  cut"  approxima- 
tions to  the  electric  field  generated  inside  a spacecraft.  (The  fields 
calculated  for  a steady  stream  of  incoming  current  are  also  valid  for 
the  case  of  an  incoming  pulse  of  current,  provided  it  can  be  assumed 
that  the  incoming  charge  distributes  itself  around  the  cylinder  very 
rapidly  in  comparison  to  the  duration  of  the  pulse.) 

Through  the  use  of  several  simplifying  parameters  and  assumptions, 
it  is  shown  that  the  problem  reduces  to  one  having  a fairly  straight- 
forward solution.  Thus,  the  primary  emphasis  of  this  thesis  will  be  a 
comparison  of  two  numerical  techniques  used  to  solve  the  differential 
equation  governing  the  problem.  These  two  techniques  are  the  use  of 
variational  calculus  and  the  use  of  Green's  functions.  Initially,  the 
differential  equation  (Laplace's  Equation)  and  the  boundary  conditions 
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are  derived.  This  equation  is  then  solved  using  separation  of  variables. 
The  results  of  this  section  are  used  as  a basis  of  comparison  with  the 
two  other  techniques.  This  is  followed  a presentation  of  the 
variational  approach.  First,  the  appropriate  minimization  functional  is 
derived.  It  is  then  converted  to  a matrix  problem  using  Simpson's  Rule 
and  a three  point  quadrature  formula.  The  matrix  equation  is  then  solved 
using  a Choleski  decomposition.  The  following  section  contains  a discus- 
sion of  the  Green's  function  method.  The  Green's  function  is  first  derived 
and  then  evaluated  numerically.  The  final  chapter  contains  a discussion 
of  the  electric  field  generated  using  representative  values  for  the  dimen- 
sions of  the  cylinder  and  the  magnitude  of  the  incoming  current. 

The  individual  chapters  contain  only  the  basic  logic  behind  the  math- 
ematical techniques  and  the  results  of  applying  them.  Detailed  derivations 
of  all  equations  can  be  found  in  the  appendices.  The  appendices  also  con- 
tain a discussion  of  the  computer  program  that  was  used  with  each  of  the 
three  mathematical  approaches. 


Physical  Analysis 


II  P'  'j.u  loj p.J  of  a Th ■ -ofot  ic a 1 Mod e l 

The  spacecraft  is  modeled  as  a non-rotating,  hollo*/,  metallic 
cylinder  with  a finite  conductivity.  If  this  cylinder  is  irradiated 
with  a uniform,  steady  stream  of  current  density,  charges  will  be 
distributed  around  the  cylinder.  This  distribution  of  charge  will 
cause  a distribution  of  potential  throughout  the  cylindrical  shell, 
including  its  inner  surface.  These  potentials,  and  the  resulting  cur- 
rent flow,  will  in  turn  induce  an  electric  field  within  the  hollow  por- 
tion of  the  cylinder.  If  this  electric  field  is  strong  enough,  undesirable 
electronic  noise  or  damage  may  result.  Knowledge  of  the  strength  of  this 
induced  electric  field  is  desirable  for  obvious  military  reasons.  The 
procedure  for  obtaining  the  electric  field  is  quite  straightforward. 

Once  the  potential  distribution  on  the  inner  surface  is  known,  the  poten- 
tial distribution  induced  in  the  hollow  portion  of  the  cylinder  can  be 
calculated.  The  gradient  of  this  distribution  will  then  give  tiie  strength 
or  the  electric  field  within  the  cylinder. 

For  a finite  cylinder  with  end  caps  the  current  flow  throughout 
the  cylindrical  shell  v/ould  have  components  both  parallel  and  perpen- 
dicular to  its  longitudinal  axis.  Also,  current  flow  through  the  end 
caps  would  have  to  be  considered.  Furthermore,  the  problem  would 
require  a solution  in  three  dimensions.  However,  by  considering  an 
infinite  cylinder,  the  end  effects  can  be  neglected,  and  all  current  flow 
will  be  perpendicular  to  the  longitudinal  axis  of  the  cylinder.  This 
simplification  reduces  the  problem  to  a manageable  two-dimensional  form. 
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The  electric  field  predicted  by  this  infinite  cylinder  model  will  approx- 
imate the  electric  field  near  the  plane  z = height/2  of  a finite  cylinder. 
The  two-dimensional  form  of  the  problem  is  illustrated  in  Fig  2-1. 


The  incoming  current  density  is  incident  perpendicularly  to  the  longi- 
tudinal axis  of  the  cylinder  over  the  interval  - tt / 2 <0 <_a / 2 . As  the  in- 
coming current  continues,  charge  will  continue  to  build  up  on  the  outside 
of  the  cylinder.  It  is  assumed  that  this  charge  lias  no  effect  on  the 
incoming  beam.  The  radially  outward  "leakage"  current  functions  to  make 
the  boundary  conditions  consistent  with  the  mathematics.  It  is  discussed 
in  the  section  explaining  the  boundary  conditions. 

Method  of  Solution 

The  current  through  an  arbitrary  surface  is 

Jf-  0-0 

S 

For  a closed  surface  this  integral  must  equal  the  rate  of  decrease  of 
charge  within  it.  With  the  charge  designated  as  q,  the  result  may  be 
written 

Jf-n  = 

s 

Through  application  of  the  divergence  theorem,  this  surface  integral 

can  be  changed  to  a volume  integral.  Also,  when  the  charge  is  represented 

as  a volume  integral  of  charge  density,  p,  Eqn  (2-2)  can  be  written 

J{v-f)jv  = -Mfdv  = 

V V 

4 


-W  « 


Since  the  surface  is  constant,  the  partial  derivative  with  respect  to 
time  can  be  brought  inside  the  integral.  Equating  integrands  yields 
the  continuity  equation 


V-T  = - ft 


3p 


(2-4) 


For  the  steady  state  condition  -p-  = 0.  For  an  isotropic  medium  (both 
regions  I and  II  are  considered  to  be  isotropic),  J = crE . Thus 


V * o*£  = o 

Replacing  E with  yields  Laplace's  Equation 

= o 


(2-5) 


(2-6) 


This  equation  is  valid  for  both  regions  of  the  cylinder.  Thus,  Laplace's 
Equation,  subject  to  appropriate  boundary  conditions,  must  first  be  solved 
for  Region  I.  This  will  yield  the  potential  on  the  inner  boundary.  Then 
Laplace's  Equation  is  solved  for  Region  II,  subject  to  the  previously  cal- 
culated values  of  (J)  at  the  inner  boundary.  This  gives  the  values  of  4> 

• t t , , 
throughout  Region  II.  The  strength  of  the  E field  in  this  region  can  then 

-V 

be  determined  from  E = -Vc|>. 

Boundary  Conditions 

Current  density,  the  electric  field,  and  potential  are  related  by 


D - cr  E - - crV  9 


(2-7) 


v<>  = - 


(2-8) 


On  the  outer  boundary,  only  the  normal  component  of  the  incoming 
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current,  E,  is  known.  The  tangential  component  cannot  be  specified 
because  the  current  flow  around  the  cylinder  is  unknown. 

On  the  inner  boundary  the  normal  component  of  7<J>  (i.e.  the  cur- 
rent) is  zero.  This  is  because  the  inside  wall  is  not  being  irradiated, 
and  all  current  flow  at  the  inner  boundary  is  tangential  to  the  bound- 
ary. There  is  no  current  flow  from  Region  I of  the  cylinder,  across 
the  inner  boundary,  into  Region  II.  As  before,  the  tangential  com- 
ponent of  V<)  cannot  be  specified  because  the  current  flow  is  unknown . 

Thus,  the  apparent  boundary  conditions  are,  for  the  outer  boundary, 


an 


a $ _ 

d~T 


O 


~ £ *0-  £ 


c & c 3jr 

2 " - a 


(2-9) 


and  for  the  inner  boundary 


OiOt  2 it 


(2-10) 


However,  the  divergence  theorem  states 


J d*  = j>  ds  <2-u) 

A S 

Since  V2 if  = 0,  it  is  obvious  that 

ji tds  • ° <2'12) 

S 

HDwever,  integrating  around  both  boundaries  with  the  ooundary  conditions 
as  specified  in  Eqns  (2-9) and  (2-10)  yields 


▼ 
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(2-13) 


(2-14) 


1PJ:  ^ O (2-13) 

<r 

Thus,  the  boundary  conditions  as  stated  in  Eqns  (2-9)  and  (2-10)  must 
be  modified  to  satisfy  the  divergence  theorem.  This  modification  con- 
sists of  adding  a radially  outward  "leakage  current"  on  to  the  outer 
boundary.  Its  value  is  taken  to  be  constant  over  the  outer  boundary, 
as  shown  in  Fig  2-1.  With  the  addition  of  this  term  the  outer  boundary 
condition  becomes 


aft 


./-J.  CO  5 -9-  + t_  j 


~~  £ -e-re 


ar 

z 


.0.  £.  3 IT* 
Z z 


(2-16) 


where  L is  a value  such  that 

bj||  d-0-  = O (2-17) 

o 

9 $ 

When  the  values  for  — are  substituted  over  the  appropriate  angular 
regions,  Eqn  (2-12)  becomes 


8 


o 


(2-13) 


3 <r 

{*}-~r  d<> 

IT 

z 


Simplifying  and  evaluating  the  integrals  yields 


-rr 

~'L  sin# 

lir 

Z 


From  which 


L = 


o 


(2-19) 


(2-20) 


Thus,  the  differential  equation  and  boundary  conditions  are, 


v2<>  * o 


(2-21) 


rs  a 


f-lo 


(cos  ■£>  - t S’ 

it  ' «*/  2 a 

_ JL  ~ c -9-  £ 3 ^r 
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I [ X So lut ion  by  Separat ion  of  Variables 


I 


Solution  in  Region  I 

Through  application  of  the  method  of  separation  of  variables,  the 
general  solution  for  Laplace's  Equation,  subject  to  the  boundary 
conditions  as  listed  in  Eqn  (2-21),  is 


(pz(r,o)  = 


X b COR  -0- 


2o-(L-«!) 


(r  * 4) 


<0 


Zi  ^ 2 / 

(-Q I h c o$  ( 2 n 9- )_ 

, hfl" o'") 
n-*l  ' ' 


/ Jn  •/*  \ 

(r  *f*) 


(3-1) 


This  equation  gives  the  potential  distribution  throughout  Region  I. 

A detailed  derivation  of  Eqn  (3-1)  can  be  found  in  Appendix  A.  Appen- 
dix B contains  the  Fourier  expansion  of  the  boundary  conditions. 


Solution  in  Region  II 

Laplace's  Equation  must  also  be  solved  for  the  potentials  through- 
out Region  II,  subject  to  Dirichlet  boundary  conditions  because  the  po- 
tentials on  the  inner  boundary  are  now  known.  The  boundary  condition 
for  Region  II  is 

•-  -f  < (3~2) 

XL 

where  f(a,0)  is  given  by  the  solution  to  Region  I at  r=a. 

The  solution  of  Laplace's  Equation,  subject  to  the  boundary  condi- 
tion as  given  in  Eqn  (3-2),  for  Region  II  is 
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(3-3) 


4jt<n*>  --  ■ VrJ— + 

2-  G ( <o  - a.  ) 

0,3  n-l  2n  + l 5 

V’  (-0  2.U r cos  (?..») 

/ . n 1>  o-  ( -v  I )1 1“"  - n''") 

US/  ' 

This  equation  is  derived  in  Appendix  A. 


'-■Field  in  Region  II 

Since  E = -Vij>,  and  the  gradient  o£  41  in  cylindrical  coordinates  is 


v?  = 


• * 

— — -o 

r 


(3-4) 


the  E-field  given  by  the  potential  distribution  in  Eqn  (3-3)  is 

ft  ^ "*•  .2"*'  2»-l 

•0  £ t)  cor, -6-  \ (-0  VI  b y ^ 

L “ “4-  cr*n-(VnM)  (£"-  aVM)  J 


. f xb(-s»v>)  jr 


rv-l  2«+'  2^-1  . H 

(-0  VJ!.  r (-emSnjcO  * 

'o-n- ( l)(b"n-  a'")  J 


(3-5) 


Appendix  C contains  a computer  program  based  on  Eqn  (3-5)  which 
calculates  the  magnitude  of  the  E-field  for  an  arbitrary  r and  G,  using 


= "VW  + La-of 


(3-6) 


Test  Solution  in  Region  I 

The  solutions  given  in  the  previous  paragraphs  represent  solutions 
for  the  problem  as  derived  in  Chapter  II.  However,  to  aid  in  the  com- 
parison of  the  results  of  the  two  other  techniques,  a simpler  outer 
boundary  condition  was  used.  This  condition  is 


3$ 

St  | 

fcj. 


- COS  ^ 


Zrr 


(3-7) 
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This  condition  has  the  advantage  of  reducing  the  Fourier  expansion  of 
f(0)  to  only  one  term. 

The  solution  to  Laplace's  Equation,  subject  to  the  above  outer 
boundary  condition  is 


4>  (<*,<>) 


coZ  -0- 


( 3-8 ) 


For  ease  of  calculation,  the  inner  and  outer  radii  will  be  2 and  4 
respectively.  Substitution  of  these  values  yields 


$(*,<►)*  •§(T  + -f)coS‘e'  <3'9 

which,  will  be  used  to  compare  the  accuracy  of  the  results  of  the  varia- 
tional and  Green's  function  approaches. 
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Solution  Procedure 


I V Variational  Calculus  Approa c 1 1 


The  principle  behind  the  variational  approach  is  that  of  minimizing 
a functional.  For  the  problem  as  outlined  in  Chapter  I the  functional, 

S , is 


s ~ 


i J ( 47  + f **do 

A O 


(4-1) 


where  the  second  integral  is  over  the  outer  boundary,  and  f(9)  is  the 
outer  boundary  condition.  There  is  no  integral  term  for  the  inner 
boundary  because  f(9)  = 0 on  the  inner  boundary.  This  functional  is 
derived  in  Appendix  D. 

Through  the  use  of  a 3-point  quadrature  formula  and  Simpson's 
Rule,  Eqn  (4-1)  can  be  written  in  matrix  form  as 

s ~ £ 5 4 ❖ — J a -f  (-»)  b (4-2) 

where  is  the  coefficient  matrix  resulting  from  application  of  Simpson's 
Rule  and  the  quadrature  formula,  and  B is  a column  vector  containing  the 
Simpson's  Rule  coefficients  for  the  mesh  points  on  the  outer  boundary. 
Taking  the  variance  of  S with  respect  to  <j>  yields 

£5  s A <*>  - 8r(-o)b  (4-3) 

Since  S is  to  be  minimized,  the  variance  must  be  zero,  so 

A <t>  - 8 <•*)  b c o (4-4) 
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Thu:;,  the  matrix  equation  is 

A $-  B -f  (&)  b 

Equation  (4-5)  can  be  written  in  the  form 

\X  $ - y 

where  L_  is  a lower  triangular  Choleski  decomposition  of  A. 
Equation  (4-6)  can  be  easily  solved  by  letting 


(4-5) 


(4-6) 


L.  x = y 


(4-7) 


and  solving  for  x using  forward  substitution.  Once  x is  known,  $ can  be 
calculated  using  backward  substitution. 

Since  the  original  problem  contains  Neumann  boundary  conditions,  it 
is  only  possible  to  specify  <j>  to  within  an  additive  constant.  This  is 
manifested  in  the  fact  that  the  coefficient  matrix,  A,  will  be  singular. 

When  A is  decomposed  via  a Choleski  decomposition,  the  lower  right  term 
in  the  main  diagonal  of  L will  be  zero.  Thus,  the  last  term  in  each  of 
the  x,  ar>d  vectors  cannot  be  solved  for  directly,  because  that  would 
entail  dividing  by  zero.  This  difficulty  can  be  alleviated  b;  only  solving 
for  the  first  N-l  terms  in  each  vector,  and  arbitrarily  setting  che  Nth 
term  equal  to  zero.  The  solution  so  calculated  will  be  correct  within  an 
additive  constant.  To  arrive  at  the  minimum  solution,  (one  which  is  ortho- 
gonal to  £ = ( 1 , 1 , 1 , . . . 1 ) , the  zero  eigenvalue-eigenvector  of  A) > it  will 
be  necessary  to  normalize  the  solution.  That  is,  total  all  the  terms,  find 
the  average,  and  subtract  this  average  from  each  term  in  the  solution  vector. 
A computer  program  for  the  variational  approach  is  listed  and  discussed  in 
Appendix  G. 


14 


Matrix  Form  of  the  Doable  Integral 

fhe  first  integral  in  F.qn.  (4-1)  can  be  converted  to  a double  integral 
ovei  r and  0,  using  the  relation  dA  — rdrdO.  The  integral  is  then  expressed 
as  a double  sum. 


r- J ~ 4'  J Jrfr/tysJ'r  d-O- 


(4-8) 


-o  r 


i j J 


(4-9) 


where  A.  and  A.  are  Simoson's  Rule  coefficients. 

1 J 

S ince 


<w)  - (i?)  + 


(4-10) 


it  will  be  necessary  to  have  difference  relations  for  both  r and  9. 
Also,  since  each  difference  relation  will  have  to  be  squared,  it  will 
be  advantageous  to  use  ones  which  don't  have  too  many  terms,  yet  are 
reasonably  accurate.  The  following  3-point  formula  are  accurate  to 
order  h2  (Ref:  8,  96). 


= ~L  (-3$  + */<$>  . - 4>a  \ 

ar  2«r  *J  ,J/ 


(4- 1 la) 


34» 

dr 


2 at 


H",i  + 


(4-1 lb) 


?)4> 

dr 


1 ( 4>  — V & . + 3 4>  \ 


2 ay 


(4-llc) 
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Eqn  (4-lla)  is  used  for  all  mesh  points  on  the  inner  boundary,  and 
Eqn  (4-lle)  is  used  for  all  mesh  points  on  the  outer  boundary.  However, 
care  must  be  used  when  approximating  the  derivative  at  the  interior  mesh 
points . 

Initially,  Eqn  (4-llb)  was  used  to  approximate  the  derivative  at  all 
interior  mesh  points.  However,  this  produced  the  scalloping  phenomena 
shown  in  Fig  4-1.  No  matter  how  fine  a mesh  was  used,  the  scalloping 
was  always  present,  and  always  over  3-point  segments.  The  reason  the 
scallops  are  present  can  be  explained  by  the  following  analysis: 

If  the  boundary  term  is  ignored,  the  functional  has  the  form  (in  one 
d imension) 

5 - if(  vf  dr 

u 


The  variance  of  S with  respect  to  $ ' is 


6 

a 


(4-13) 


If  S is  a minimum,  6 S = 0.  Also,  if  the  curve  (a,b)  is  approximated  by 
two  piecewise  continuous  functions  in  the  intervals,  (a,c  ) and  (c  ,b), 
Eqn  (4-13)  can  be  written 


(4-14) 
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IE  both  integrals  are  integrated  by  parts,  this  equation  bt comes 


C b 


dr  - O 


(4-15) 


Evaluating  the  Cirst  two  terms  and  combining  the  integrals  yields 


<(<>,<  - + / - 

b 

//  Cxj>"  = O ^ 


(4-16) 


Due  to  the  boundary  conditions,  the  terms  evaluated  at  "a"  and  "b"  are 
both  zero.  Also,  since  the  function  in  the  interval  from  a to  c to  b is 
continuous,  Eqn  (4-16)  can  be  written 


(4-17) 


Thus,  since  continuity  of  the  first  derivative  is  not  enforced  across  each 
segment,  the  functional  has  been  given  the  freedom  to  absorb  some  of  the 
"cost"  of  minimization  by  letting  the  slope  be  discontinuous  across  each 
piecewise  continuous  segment.  The  reason  the  scallops  occurred  in  groups 
of  three  points  is  because  both  Simpson's  Rule  and  a 3-point  quadrature 
formula  were  used  in  the  numerical  approximation,  implying  a quadratic 
function  over  a 3-point  range. 


To  eliminate  the  scalloping,  Eqns  (4-lla)  and  (4-llc)  were  used  to 
approximate  both  the  left  and  right  derivatives  at  the  mesh  points  which 
are  junctions  for  each  of  the  3-point  segments  (e.g. , mesh  point 


i.J 


in 


Fig  E-l).  This  procedure  more  accurately  calculates  the  actual  value  of 
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tho  derivative  on  either  side  of  the  break-point,  so  that  large  slopes 
(scalloping)  will  be  suppressed. 


Eqn  (4-llb)  can  be  used  to  approximate  the  derivative  at  the  middle 
mesh  point  of  each  3-point  segment.  In  this  case  continuity  is  enforced 
through  the  use  of  Simpson's  Rule.  Use  of  the  approximation  formulae 
in  this  manner  produced  a smooth  curve  for  the  approximation  in  the  radial 
direction. 

Similarly,  care  must  be  used  when  approximating  the  angular  deriva- 
tives of  Eqn  (4-10).  Initially,  angular  derivatives  at  all  mesh  points 
were  approximated  with  Eqn  (4-llb).  As  can  be  seen  in  Fig  4-2,  the  results 
of  this  approximation  bracket,  but  never  approach,  the  analytic  answer. 

The  values  along  the  even  numbered  radials  are  always  greater,  and  the 
values  along  the  odd  numbered  radials  are  always  less  than  the  analytic 
answer.  Again,  this  results  from  an  improper  approximation  of  the  first 
derivative  in  the  angular  direction.  This  was  corrected  by  using  Eqns  (4-lla) 
and  (4-llb)  to  approximate  the  left  and  right  angular  derivatives  at  all  mesh 
points.  Again,  an  additional  factor  of  \ was  included  in  both  Eqns  (4-lla) 
and  (4-llb)  to  account  for  the  use  of  both  left  and  right  derivatives.  Use 
of  the  approximation  formulae  in  this  manner  produced  a smooth  curve  for  the 

•r 

approximation  in  the  angular  direction. 

Ma t rix  Form  of  the  Single  Integral 

The  second  integral  in  Eqn  (4-1)  can  be  converted  into  the  following 
summation: 


J-2  ■fi&J'rcl-G-  = — £ 8.  U-o-). 
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1 


i 


No  d L £ ference  formulae  are  needed  because  there  are  no  derivatives  to 
approximate.  In  fact,  the  only  non-zero  terms  in  the  left  hand  side 
vector  will  be  those  terms  corresponding  to  a value  of  •>  on  the  outer 
boundary . 


Summary 

Initially,  the  boundary  conditions  were  enforced  with  Lagrange 
Multipliers,  without  considering  the  niturai  boundary  conditions  of  the 
functional.  This  produced  a scalloped  curve  with  endpoints  at  or  near 
zero.  (This  initial  attempt  is  discussed  in  Appendix  F.)  When  the 
natural  boundary  conditions  were  included  in  the  functional,  the  curve 
was  "pulled  up"  so  that  its  endpoints  more  closely  approximated  the 
analytic  values.  However,  the  scallops  were  still  present  no  matter 
how  fine  a mesh  was  used.  These  scallops  were  present  because  continuity 
of  the  first  derivative  was  not  rigidly  enforced  in  the  numerical  approx- 
imation. When  both  left  and  right  radial  derivatives  of  the  mesh  points 
which  were  junctions  of  the  3-point  segment  were  used,  the  numerical 
solution  did  indeed  become  a smooth  curve.  However,  it  still  did  not 
satisfactorily  approximate  the  analytic  solution.  Finally,  both  left 
and  right  derivatives  in  the  angular  direction  were  approximated  for  all 


mesh  points.  iCs  illustrated  in  Fig.  4-3,  the  numerical  solution  approaches 
the  analytic  result  as  the  mesh  becomes  finer.  However,  for  a 9 x 32 
mesh  it  was  necessary  to  use  168  K of  computer  storage. 
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V So  Lut  Ion  I’. sing  Green 'a  Func  t i o n s 


In  this  chapter  the  equations  developed  in  Chapter  1 are  solved 
using  a Green's  function.  First,  a solution  for  if>  is  derived  in  terms 
of  a Green's  function.  Second,  an  iterative  scheme  for  obtaining  the 
values  of  $ is  presented. 


Reduction  to  Integral  Equation 
The  basic  equation  is 


with  boundary  conditions 


9 <y>  s O 


s4> 

or 


= O 


r.-a 


hrj 

1*9  0 


- Y ("&) 


The  Green's  function  expression  for  Eqn  (5-1)  is 


(5-1) 


V Gi(3 ?,F)  - /(*-*') 


(5-2) 


where  x is  the  observer's  point  and  x'  is  the  source  point,  as  measured 
from  the  origin.  Since  the  Green's  function  only  depends  on  the  difference 

between  the  two  vectors,  x and  x',  it  can  be  written  as  G(x-x'). 

Multiplying  Eqn  (5-1)  by  G(x  - x')  and  Eqn  (5-2)  by  4>(x)  yields 

r*  j \ i j r* 


G(x-  k')v*t>(x)  = o 


(5-3) 


and 


yzo>(x-x')  ~ <f(3T->r)  4>(*) 


(5-4) 
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Subtracting  Equ  (5-J)  from  F.qn  (5-4),  and  integrating  over  the 
volume  of  the  cylinder  results  in 

/[-&{?- S')  + «(.?)V‘£i;A;-3?')]dV  = <*>(*')  (5-5) 

V 

If  one  of  the  "del"  terms  is  factored  out  of  the  brackets,  the  volume  integral 
can  be  converted  to  a surface  integral  through  application  of  the  divergence 


theorem. 


*')  - v4>(^)]#  A*  - <>(*')  (5_ 


6) 


Since  the  Green's  function  is  symmetric,  the  observer's  point,  x,  and  the 
source  point,  x',  can  be  interchanged.  This  change  results  in 


<Hx)  VG|(X~*)  - 6 {*'-*)  ($')]»  cl  a 


If  a new  function,  i'(x),  is  defined  as 


Vex)  = -Je, v'<M?')-da' 


Eqn  (5-7)  can  be  written  as 


4><?)  = f it)  + J <H*‘)  v'c,  (T-  X ) • da' 


Equations  (5-8)  and  (5-9)  can  be  solved  by  iteration,  viz., 


4>0(*>  - 7 (x) 


4>n  (*)  - U/(K)  + J (*)  V& (*-*)•  da' 


(5-7) 


(5-8) 


(5-9) 


(5-10a) 


(5-10b) 
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where 


(5-11) 


I?- 


and 


vc,(r-7) 


/ —A 

X - * 

* \x'-*\z 


(5-12) 


Eqns  (5-11)  and  (5-12)  are  derived  in  Appendix  II.  Also,  a more  rigorous 
derivation  of  the  integral  equation,  Eqn  (5-7)  can  be  found  in  Appendix  I. 

The  iteration  scheme  is  explained  in  the  following  six  steps. 

1.  The  cylindrical  shell  is  divided  into  an  even  number  of  equally 
spaced  angular  intervals.  Since  Eqns  (5-8)  and  (5-9)  only  contain  surface 
integrals,  there  is  no  need  for  interior  mesh  points.  (This  is  one  advan- 
tage of  using  Green's  functions  - the  number  of  mesh  points  can  be  dras- 
cicaLly  reduced,  thus  saving  computer  storage  space.) 

2.  Equation  (5-8)  is  used  to  calculate  the  value  of  tj/(x)  for  each 
mesh  point  on  the  inner  boundary  due  to  contributions  from  every  mesh 
point  on  the  outer  boundary. 

3.  Equation  (5-8)  is  used  to  calculate  the  value  cf  ii(x)  for  each 
mesh  point  on  the  outer  boundary  due  to  contributions  from  every  mesh 
point  on  the  outer  boundary. 

Note:  In  both  steps  2 and  3 there  is  no  contribution  from  the  inner 

boundary  because  V({>  * da  = 3i}>/3r  = 0 on  the  inner  boundary. 

4.  Equation  (5-9)  is  used  to  calculate  the  value  of  4>  for  every  mesh 
point  on  the  inner  boundary  due  to  contributions  from  all  mesh  points  on 
both  boundaries.  The  correspond ing  value  of  \ !>  is  then  added  to  the  sum  of 
these  contributions. 

5.  Step  4 is  repeated  for  all  mesh  points  on  the  outer  boundary. 
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6.  Steps  4 and  5 are  repeated  until  all  values  of  tj>  converge  to  within 
some  tolerance  limit. 

A computer  program  implementing  the  above  iteration  scheme  can  be  found 
in  Appendix  J. 

The  values  of  tp  calculated  by  this  method  are  in  agreement  to  three 
significant  figures  with  the  values  predicted  by  the  test  solution  of  Eqn 
(3-9).  In  this  calculation,  128  mesh  points  were  used  around  each  boundary 
and  the  values  of  the  conductivity  and  incoming  current  were  set  to  unity. 
The  computer  program  was  then  modified  to  accommodate  the  boundary  condi- 
tions as  listed  in  F.c;n  (2-21).  Once  again,  the  values  of  d>  were  in  agree- 
ment to  three  significant  figures.  A comparison  of  the  results  was  not 
plotted  because  the  curves  would  be  indistinguishable. 


Solution  in  Region  II 

The  procedure  is  the  same  as  that  used  in  the  previous  section  up  to 
Eqn  (5-7).  On  the  inner  boundary  Vj>(x')  = 0,  so  Eqn  (5-7)  can  be  written 


-•  / V&cx'-x)  <£(«')•  da' 


(5-13) 


where  VG(x'  - x)  is  given  in  Eqn  (5-12)  and  the  values  of  <Kx')  have  been 
calculated  by  the  method  presented  in  the  previous  section. 

Thus,  it  is  not  necessary  to  iterate  a solution,  just  sum  the  contri- 
butions from  all  points  on  the  inner  boundary.  The  summation  expression 
for  Eqn  (5-13)  is 


N 

j COS  dn  <K#), 

*<*  > - — L 

n*|  1 n 


(5-14) 


where  y is  the  angle  between  x -x  and  the  outward  unit  normal  at  x . 


1 


A 
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The  potentials  obtained  with  this  formula  are  in  agreement  to  two 


significant  figures  with  the  values  predicted  by  Eqn  (3-3). 

Calculatio n of  the  11-Field 

It  is  not  necessary  to  calculate  the  values  of  <1  throughout 
region  II.  Eqn  (5-13)  can  be  used  to  calculate  the  electric  field 
directly.  Since  E = — V4> » Eqn  (5—13)  can  be  written  (after  the  dot 
product  operation  is  performed) 

— <!>(*)  “ -vf  ^ <£>(£')  ds  (5-15) 

5 

The  gradient  operation  is  with  respect  to  x,  not  x',  so  Eqn  (5-15)  becomes 


-Ji?C 

S 


r® 

/<)& 1 \ * , 

/36)  £ 

[3t 

(^i  j T 

V* 

ds 


(5-16) 


Summary 

The  Green's  function  method  approximated  the  values  of  <1  as  predicted 
by  the  analytical  formulas  of  Chapter  III  to  at  least>  two  significant 
figures.  Thus,  the  Green's  function  method  gave  an  independent  verifica- 
tion of  Eqns  (3-1)  and  (3-3).  Time  did  not  permit  a computer  implementa- 
tion of  F.qn  (5-16). 
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VI  ' P r ' ' s 0 n r -i  *-  i v e V 1 u c s o t~  the 

, -> 

.Magni  tude  of  the  E-Field 
- 


In  this  chapter  the  electric  field  induced  inside  the  hollow  portion 
of  the  cylinder  is  calculated  using  Eqns  (3-5)  and  (3-6).  The  following 
representative  values  are  used  for  the  input  parameters: 

Outer  radius:  .5  meters 

Cylinder  thickness:  varied  from  .254  to  2.54  millimeters  (10  to  100  mils) 

Conductivity  (Aluminum):  3.82  x 107  mhos/meter 

Incoming  current  density:  .01  amps/meter2 

Fig  6-1  illustrates  how  the  maximum  strength  of  the  electric  field 
varies  with  the  thickness  of  the  cylinder. 

Fig  6-2  illustrates  the  magnitude  of  the  electric  field  as  it  varies 
with  angle.  Each  curve  plots  the  field  along  an  equi-radial  arc.  Since 
the  electric  field  is  symmetric  with  respect  to  the  line  0 = 0,  only  the 
fields  in  the  upper  half  of  the  cylinder  are  plotted.  As  can  be  seen 

• • “V  . 

from  the  figure,  the  maximum  value  of  the  E-Field  occurs  at  the  inner 
boundary  when  0=0. 

Both  figures  indicate  that  for  an  incoming  current  density  of  .01 
amps/meter2,  the  strength  of  the  electric  field  inside  the  cylinder  is  on 
the  order  of  10-7  volts/metcr.  Eqn  (3-5)  indicates  that  the  magnitude  of 
the  electric  field  varies  directly  as  the  incoming  current.  Thus,  to  start 
getting  appreciable  fields  inside  the  cylinder,  the  incoming  current  density 
would  have  to  be  on  the  order  of  10^-10°  amps/meter2. 
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Fig  6-1.  Plot  of  Ik  I vs  Cylinder  Thickness 

'max 

for  an  Aluminum  Cylinder  with  a Radius 
of  .5m. 
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A comparison  of  the  variational  computer  program  with  the  Green's 
function  program  shows  that  the  variational  program  requires  less  execu- 
tion time,  whereas  the  Green's  function  program  requires  less  computer 
storage.  Thus,  for  a problem  which  cannot  be  solved  analytically,  the 
best  numerical  approach  would  depend  on  whether  computer  storage  or  execu- 
time  is  more  critical  to  the  user. 

Lastly,  given  the  representative  values  for  the  cylinder's  dimensions 
and  the  incoming  current  density  as  listed  in  Chapter  VI,  the  magnitude  of 
the  induced  E-Field  is  on  the  order  of  10-7  volts/mcter. 

Recommendations 

1.  The  computer  programs  developed  in  Chapters  IV  and  V only  calculate 
the  potential  distribution.  It  would  be  interesting  to  develop  them  one 
step  further  to  calculate  the  E-Field , and  then  Compare  the  two  numerical 
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methods . 


2.  To  improve  the  prediction  of  the  electric  field  inside  a space- 
craft, the  problem  should  be  solved  for  a finite  cylinder  with  endcaps. 
One  possible  model  is  a thin  ellipsoid.  This  would  alleviate  the  problem 
of  calculating  the  current  flow  across  the  edge  of  the  endcaps. 
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Appendix  A 


12ll£  ' v:l  *■  1 ou  _?  1 the  Solution  to  V 2 0 = 0 
t's  i ng  S >'p,n _n  t io n o f Variables 


Ti\e  method  of  separation  of  variables  is  used  to  solve  V2<j>  = 
The  boundary  conditions  are  then  expanded  in  a Fourier  series  to 
an  infinite  sum  expression  for  {> ( r , 0 ) . 


So lu tio n in  Region  I 


= o 


Assume  - Rfr)  O(0).  In  cylindrical  coordinates 

v?24>  - f±  2.  (r  JL  \ s JL  1 

[riirl  / ' V a <>> 

After  expanding  and  multiplying  by  r2, 

d /rJA)  ~ _ 

Jr  V dr  / 

Solution  if  A2  = 0 


Rd 


r d 

R Jr 


I s ^ 


e 


Solving  for  R: 


(r  0 

civ  ^ /-I  v* 


cjr  / 


r - A 

r dT 


F>  ~ A,  In  f + B 
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0. 

obtain 


(A-l) 

(A-2) 

(A- 3) 


(A-4) 

( A— 5 ) 

( A-  6 ) 


Solving  for  0: 


do* 


(A-7) 


e = c0o  * dq 


(A-8) 


* (*0  U r + 30)  ( C.+  + D.) 


(A-9) 


Since  4>  is  periodic,  C0  must  be  zero.  Absorbing  Dc  in  the  remaining 
constants  yields 


<P  = In  r -4*  S4 


(A-10) 


Solution  if  A?  is  Positive 


Solving  for  R: 


rf;(^  ~ = 


( A— 11) 


This  is  Euler's  Equation,  which  can  be  solved  by  letting 


r - e 


(A-12) 


ciR  a 


dr  dz  cir 


(A-13) 


The  solution  is 


R a A r*  + B r' 


( A- 1 4 ) 


Solving  for  0: 
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o 


( A-l  5 ) 


.It  _ 

c[0 

k\0Z 


-{-  A©  = 


(~)  ~ C Coil  A*G-  + O Sin  nO- 


(A- 16) 


For  solutions  periodic  in  2rr,  X = n,  where  n = 1,2,3,...  . Thus  the  solu- 
tion for  a positive  A2  is 


V - + 0nY'V'j(C?1 04^  -h  Dn3i n n <$j  h»>,2,— 

The  solution  when  A2  is  negative  can  be  ignored  because  it  contains 
sinh  6 and  cosh  0 terms,  which  are  not  periodic. 

Thus,  the  general  solution  is, 

4>  = Aa  I*  v + 5to  + 


( A- 1 7 ) 


H [(^n'*n  + BnV  "jcos  + (Cn  rn  -1-  Dn  r'n}  S,in  r.^]  (A  18) 

nsi 

The  boundary  conditions  are 


ar 

£>t 


-V?-)  ~ o 


= -f  (-0.)  = < 


-_I_ 

<T  TT 


_TT  . TT 

— C.-0-  C 

a ~ ~2 


H A -o-x.  3jr 

2 " a 


Applying  the  first  boundary  condition  yields 


4 a 

r\ 


A0=  o 


r>- 1 -n-| 

- hn*  ~ o 


Cn“" 


■I  _ -n-l 

- D^a  =r  O 


(A-19a) 

(A- 19b) 
(A- 19c ) 
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Applying  the  second  boundary  condition  yields 

cO  ft* 

^ (n/tAb  -nB„b  jcos  n#  -f  (nC^b  ~ nDnb  Js'irtn-i}'  -f(O) 
n = i 


(A-20) 


The  Fourier  expansion  Cor  f(0)  is  given  in  Eqn  (B-9),  Appendix  B.  Since 
there  are  no  sine  terms  in  the  expansion 

n-i  -n-> 

•AC.b  - nDnb  = 0 


( A — 2 1 ) 


Solving  this  equation  simultaneously  with  Eqn  (A-19c)  yields 

Cn  — Dn  - O 

Equating  coefficients  of  the  cosine  terms  in  Eqns  (A-20)  and  (jj— 9)  for 
n = 1 yields 


(A-22) 


A ~ B,  b 


-2 


(A-23) 


2 <r 


Solving  this  equation  simultaneously  with  Eqn  (A-19b)  for  n - 1 yields 

A X b (A-24a) 

1 * a2) 


8.  s 


(A- 24b) 


Z<r(g-c?) 

Equating  coefficients  of  the  cosine  terms  in  Eqns  (A-20)  and  (B— 9 ) for 
n = 2,4,6,  .. . yields 


n-2 

_ (-1)  21 


r\Anb  - nB.b 


n.  z.v.fc. 


(A-25) 


<rTr  (r^-|) 

Solving  this  equation  simultaneously  with  Eqn  (A-19b)  for  n>l  yields 
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(-0  2i  r( 

ncr-r,*^-/) 'J- 

/ |\  & pT  2»\ 

* ; .<  J_  c.  o 

no- -rr  (•rr-i)(^n-  ?"[ 


tt3  2,  V,  6, 


(A-26a) 


n--  2 ,'/,&,••• 


(A-26b) 


The  coefficients  of  cos(nO)  are  zero  for  n = 3,5,7,...  . Thus  An=Ba=0 
when  n = 3 , 5 , 7 , . . . . 

The  equation  for  i)(r,0)  can  now  be  found  by  substituting  Eqns  (A-19a), 
( A-24) , and  (A-26)  into  Eqn  ( A — 1 S ) : 


4>  <**,-©■) 


Zb  CCS  -fl-  / y.  * „2  V 
2 '(&-«’)  ' > / 


M)'a  21b  cos_n«-  /f"+  al"  \ 
n <r  li-  ( v A- 1 ) ( £ a‘ ')  \ ' r"  ) 


«=2,v,6,—  (A~27) 


Equation  (2-27)  can  be  rewritten  with  an  infinite  sum  by  substituting 


2n  for  n. 


<b  (V,  £M  r.  — ZJ (V  + .&!  \ 


ft-1  ?.r\  + l 


V (-0  X b cos(gtvi»  /r2n  -V  AVw  \ 
Z.  n<r ir  7}  ( Z")  ' / 

n=i  * ' x ' 


( A-28) 


Equation  (A-28)  gives  the  potentials  throughout  Region  I.  The  constant 
Bq  in  Eqn  (A-18)  has  been  set  to  zero.  Since  V <f>  is  the  torn  of  interest, 
any  constant  term  would  vanish  when  the  gradient  is  calculated. 


Solution  in  Region  II 


The  equation  to  be  solved  is 
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= o 

with  boundary  condition  .(>  (a,0)  - f(a,0)  where  f(a,0)  is  given  by  the 
solution  to  Region  L. 

Using  the  same  procedure  as  before,  Clio  general  solution  is 


(A-29) 


4*  = A0  In  t + 0o  •*•  ) f/Atir"n+'  8nr  )c«s  t*  n-oj  ^ 

its  I 

Since  f>  must  be  bounded  at  r = 0,  Ao~Bn~Dn=0,  A1-s°>  the  additive  con- 
stant, B0,  is  set  to  0.  Since  the  gradient  of  (J>  is  desired,  BQ  would 
drop  out  anyway.  Thus,  the  general  solution  is  reduced  to 


(A- 30) 


r°C.o5  + Cn  V*n  S»n 

rul 

Applying  the  boundary  condition  at  r = a yields 

oo 

^ co3(,v0)  + C.n«n sin  (n£-)j  = 

n.s  i 

2 J5. L ,«-•  Jn  + I 

c»sO  , V (-0  J !o  go5(2w-9)  (2«  n) 

« j V‘  4 ” “n  cr  -- r ( ^n*.  /)  ( ^^7“ 


Since  there  are  r.o  sine  terms  in  the  right-hand  side,  Cn  = 0. 
Equating  the  coefficients  of  cos  0 for  n = 1, 

rt  ..  _ I V>  CX 

< or(^’-  a2) 


which  simplifies  to 


A = _li_ 

' o(b2-  a!) 

Equating  the  coefficients  of  cos  0 for  n = 2,4,6,..., 


(A-31) 


(A-32) 


( A- 3 3 ) 


(A-34) 
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(A-35) 


cO 

CO 

I 

-E 

nel 

he  t 

Solving 

for  the  constant 

n-l  Jn+I  - „ 

(-D  21 '»  *’ 


;v”) 


# 2*+l 

(- ! ) 21  b 


n cr  fl-  »)  ( b*".  aVn) 

When  n — 3,5,7,...  the  constant  coefficients  are  zero. 

Substituting  Eqns  (A-34)  and  (A-36)  into  Eqn  (A-31)  yields  the 
potential  distribution  for  Region  II, 


(A-36) 


» / 

c'  ( b?'- 

cO 

n-l  2n  + l 

V . 

C-i)  21  b 

z_ 

n-t 

+ 


CCS  (2  r\&) 


J 


(A-37) 
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Appendix  B 


Four iii  i pans  ion  of  the  Boundary  Conditions 


In  Appendix  A a •■••lution  for  V^<j>  — 0 was  derived  in  terms  of  an 
infinite  sum  of  situ 1 cosine  terms.  The  boundary  conditions,  repre- 

sented as  f ( 0 ) , wer<  ' n expanded  in  a Fourier  series  to  determine  the 
coefficients  for  the  -me  and  cosine  terms  in  the  solution.  The  purpose 
of  this  appendix  is  (-•  present  the  Fourier  analysis  used  to  determine 
those  coefficients. 

The  solution  at  > ^ b can  be  written  in  the  form,  (Ref:  i;qn  (A-20)), 

oO 

y [Enc-Cn*>  + Fn  sin  C«^)]  C^-)  (B_l) 

n*  I 

where 


•f(o) 


f l + - *) 
< 


O 


I. 

•r 


"2T  f1  i 2 
z ~ “ z 


2 “ ~ Z 


Since  f(9)  contains  ee1/  even  functions,  its  Fourier  expansion  will  not 
have  any  sine  terms.  nonce  Fn  = 0 in  Eqn  (B-l). 

The  expression  ls 


S.  J J (■&)  COB  JfO-  he  ” 


( B-2 ) 


where  2L  = Period.  F"*  this  problem  the  period  is  2ir . 

When  the  equations  f"'  t(®)  ere  substituted  in  Eqn  (B-2)  and  integrated 
over  the  appropriate  I units,  the  following  equation  results, 
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A 


When  Eqns  ( B — 1 ) and  B-S)  are  combined  and  simplified,  the  Fourier 
expansion  for  f(0)  over  the  interval  0 to  2n  can  be  written 

£.t&\  - l££i±  + & n=  2,*v,  6,-  (b-9) 

' Zar  5 o” it  (n4-  J) 


This  expansion  was  used  with  Eqn  (A-18)  to  evaluate  the  constants  in  the 
general  solution  for  <i>(r,9). 


This  program  computes  the  magnitude  of  the  electric  field  at 

various  points  throughout  the  hollow  portion  of  the  cylinder.  It  is 

essentially  a program  of  Eqn  (3-5). 

Cards  1 - 28:  Self  explanatory. 

Card  29:  The  value  of  the  inner  radius  is  calculated. 

Cards  33,  34:  These  cards  convert  the  angular  and  radial  increments 
into  DO  LOOP  indices. 

Card  35:  A DO  LOOP  is  initialised  which  spans  all  the  angular  increments. 

Card  36:  THETA  is  the  angular  displacement  of  each  radial  line. 

Card  37:  A DO  LOOP  is  initialized  which  spans  all  the  points  along  a 
given  radial  line. 

Card  38:  R is  the  radial  displacement  of  the  point  being  calculated. 

Card  39:  This  card  ensures  that  all  points  lie  in  the  hollow  portion 
of  the  cylinder. 

Card  40:  THETA  is  converted  into  radians. 

Cards  41  - 49:  The  radial  component  of  the  electric  field  is  calculated 
as  given  in  Eqn  (3-5).  The  summation  is  continued  until  two  successive 
values  differ  by  less  than  .0001. 

Cards  50  - 58:  The  angular  component  of  the  electric  field  is  calculated 
as  given  in  Eqn  (3-5).  The  summation  is  continued  until  two  successive 
values  differ  by  less  than  .0001. 
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Card  59:  The  magnitude  of  the  E-Fie.ld  is  calculated. 


Card  60:  THETA  is  converted  back  to  degrees. 
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Appendix  D 


Derivation  of  the  Minimization  Functional 


This  appendix  contains  the  derivation  of  the  functional,  Eqn  (4-1), 
which  was  minimized  in  the  variational  approach.  The  discussion  method 
used  will  be  to  first  present  a general  functional  for  the  problem  being 
solved.  This  will  be  followed  by  a presentation  of  a specific  functional 
compatible  with  the  given  boundary  conditions  as  listed  in  Chapter  II. 

It  will  then  be  shown  that  this  functional,  J [o]  , is  a minimum  by  adding 
another  non-zero  function,  v,  to  * , and  showing  that  J [i+v]  > JW  . 

The  general  problem  is  listed  in  Eqn  (D-l).  The  equation  is  inhomo- 
geneous and  must  satisfy  mixed  boundary  conditions. 

v*4>  * s <D"U 


On  the  Inner  Surface, 


<T(  cj) 


On  the  Outer  Surface, 


3* 


ar 

r=b 


+ 


The  applicable  functional  is 


3* 


T[40  « j/((V<i>/-2<tS]cU  + 

A 

*“  24>s.)^3  + 2)ds 
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For  the  problem  as  derived  in  Chapter  II,  S = Oj  = = gj  = 0>  and 

; 2 = t( 0 ) • So  Eqn.  (0-2)  becomes 


I tw  = zf(v<b)*<iA  + ~j-z< y;«»bdo 


(D-3) 


where  the  second  integral  has  been  converted  to  an  integral  over  0 by 
use  of  the  relation  ds  = rd8=bdQ. 

It  will  now  be  shown  that  if  a non-zero  function,  v,  is  added  to  4> , 
then  jOH  >'^  D1  ‘ (both  ij>  and  v belong  to  the  class  o(  functions  which 
has  a continuous  second  derivative.)  If  this  is  true,  then  F.qn.  (D-3)  is 
the  correct  functional  to  be  minimized. 


J = ~ j [ (v4>)2  4 2 73  • w 4 (wfjdA  4 


’ff-zW+Wfwbc le- 


( D-4 ) 


= 4 jrfiwfdA  4/7-3*  w dA  - 


(D-5) 


The  second  term  is  positive  for  any  non-zero  v.  Thus,  if  it  can  be  shown 
that  terms  three  and  four  are  zero  for  any  v,  then  J [$+vJ  is  greater  than 

jW- 
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An  identity  from  vector  analysis  gives  the  relation 


_»  _s  -s 

V-(W  A)  ~ w V-  A 4-  A * Vvv 

By  lotting  v = w and  V <j)  = A , the  above  equation  becomes 

v-  (vv4>)  = + V4>  • V V 

Thus 

V 4> • V V - V • ( vv40  - '/  V-  V 4> 


Integrating  Eqn  (D-8)  over  the  area  of  interest  yields 

j ( V V V/  cl  .1  - J 7 ' (V  V 4>)  cl  A ~ J ( V v * v9)  cl  A 1 

/I  <4  /A 

Through  application  of  the  divergence  theorem,  the  first  integral  on 


the  right  side  can  be  expressed  as 

fs7-(vv<P)d  A-JVfgJfl  +/vf£*  <»-“> 

A ^inner  °oMie< 

Substituting  back  into  Eqn  (D-9)  yields 

/(vf> -VV)c|A  -/v||c!s  +J(wv<t>)dA  (D-U) 

A A 

Note:  = 0 on  the  inner  boundary. 

c)  “ d 1 

As  stated  in  the  discussion  after  Eqn  (D-5),  it  must  be  shown  that 


J(v<$  •w)dA  - Jv-f(-O)  b d-fr  = O 

A -0- 
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Substituting  Eqn  ( D— 11)  into  first  integral  yields 


Jv~~bd0-  — JVv*  V4>  d4  bd#  = O (D_l3) 

-£>  ■& 

where  the  integral  over  S has  been  converted  to  an  integral  over  6 using 
ds  = rdO  = bdO . 

Equation  (D-13)  can  be  rewritten  as 

jv  d 4 + Jv  [ - -f  (o) 

A 4- 

The  first  integral  is  zero  from  the  initial  equation  V2iJ  = 0.  The  second 
term  is  zero  from  the  outer  boundary  condition.  Thus  J ^f+vj  JptTj  for  any 
non-zero  v.  Hence  Eqn  ( D—  3 ) is  the  proper  minimization  functional  for  the 
given  equation  with  Neumann  boundary  conditions. 


bd#  = 0 


1 

Append ix  E 

Reduction  of  the  Coefficient  Matrix 


This  appendix  contains  an  explanation  of  the  use  of  symmetry  in  reducing 
the  size  of  the  variational  matrix  problem.  Also  included,  is  a discussion 
of  an  attempt  to  further  reduce  the  matrix  through  "decoupling". 

Use  of  Symmetry 

The  coefficient  matrix,  A,  generated  by  Eqns  (4-9),  (4-10),  and  (4-11) 
becomes  very  large  as  the  mesh  becomes  finer  and  finer.  However,  it  is 
possible  to  reduce  the  size  of  this  matrix  by  taking  advantage  of  the  symmetry 
of  the  problem. 

Figure  K-l  illustrates  a 5 x 8 mesh  over  the  shell  portion  of  the  cylinder. 
A mesh  this  size  was  chosen  for  explanation  because  it  illustrates  the  salient 
features  of  the  matrix  reduction  without  becoming  too  cumbersome.  The  points 
of  discussion  can  be  generalized  to  any  size,  matrix  which  has  an  odd  number 
of  radial  mesh  points,  NR,  and  an  even  number  of  angular  mesh  points,  NA.  The 
matrix  generated  by  Eqns  (4-9),  (4-10),  and  (4-11)  has  the  form 


d fgooogf 
fefhoooh 
gfdfgooo 
oh  fefhoo 
oogfdfgo 
ooohfefh 
gooogfdf 
fhooohfe 


(E-l) 


The  terms  d,  e,  f,  g,  and  h are  block  submatrices  of  size  NR  x NR. 

Since  the  incoming  current  density  is  the  same  above  and  below  the  line 
9=0,  the  values  of  the  potential  will  be  the  same  for  those  mesh  points 
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Fig  E-l.  Sample  Mesh  for  Calculating  (J>  in  the  Cylindrical  Shell. 


which  are  symmetric  about  the  line  0 = 0,  (i.e.,  <j> . . . = A . . ).  Thus, 

hj  + l i.J-1 

the  number  of  mesh  points  actually  calculated  can  be  cut  approximately  in  half. 
The  matrix  form  of  the  problem  is 


A new  variable  f,  is  defined  a£ 


A <i>  - 6 


V = 6 4> 


( E— 2 ) 


(E-3) 


where 


*•«  = i <Ai  + \->)  K'~  j= 


V - cb  - cj>.  . 


Thus  £ ~ fl  T , and  Eqn  (E-2)  becomes 


j = K-  -A 


AOW  = S 


Pre-mu It  ip  lying  Eqn  (E-4)  by  Q yields  a new  matrix  problem 


(E-4) 


The  form  of  0 is  such  that 


OA  = OB 


O S = 


( E— 5 ) 


(E-6) 


That  is,  all  the  terms  of  _B  lying  below  the  axis  of  symmetry  become  zero.  As 
can  be  seen  from  Fig  E-l,  the  values  of  B which  are  equal  are  B„,  = Bg,  B^  = B?, 
and  B^  = Bg . Since  Bg,  By,  and  Bg  all  lie  below  the  axis  of  symmetry,  they 
must  be  set  to  zero.  Thus,  it  is  necessary  to  find  an  Q 1 matrix  which, 
when  operating  on  J5,  sets  Bg  = B?=  Bg=  0. 
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II 


An  0 * matrix  which  satisfies  this  condition  is 


r- 

— 

r- 

r-  "i 

i 

0 

0 

0 

0 

0 

0 

0 

B1 

B? 

B? 

0 

~2 

0 

0 

0 

0 

0 

*-• 

0 

0 

0 

0 

0 

0 

0 

0 

0 

'2 

0 

0 

0 

«4 

0 

0 

0 

0 

1 

0 

0 

0 

Bs 

(J 

0 

0 

1 

0 

-1 

0 

0 

■V 

0 

0 

0 

1 

0 

0 

0 

-1 

0 

»7 

B 

0 

0 

1 

0 

0 

0 

0 

0 

-1 

0 

L. 

— 

8 

( E— 7 ) 


_1 

Of  course  the  0 matrix  can  be  expanded  to  accommodate  any  B vector  gen- 
erated by  an  even  number  of  angular  mesh  points. 

The  effect  of  the  operations  0 A 0 and  03  can  be  calculated  once  0 1 is 
known.  The  matrix  generated  is  significantly  smaller  than  the  original  coef- 
ficient matrix. 


where  k = e+h- 


d 

2f 

2g 

0 

0 

0 

0 

0 

2f 

2k 

2 f 

2h 

0 

0 

0 

0 

2g 

2f 

2d 

2f 

2g 

0 

0 

0 

0 A 0 = 

0 

2h 

2f 

2k 

2f 

0 

0 

0 

0 

0 

2g 

2f 

d 

0 

0 

0 

0 

0 

0 

0 

0 

X 

X 

X 

0 

0 

0 

0 

0 

X 

X 

X 

! 

0 
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0 

0 

0 

X 

X 

X 

submatrices , 

ind icated 

with 

an 

are  \ 

(E-8) 


the  corresponding  terms  in  B will  be  zero  through  the  operation  0 B. 


OB  = 


2B„ 

2B' 

2B.; 


(E-9) 


Thus,  the  matrix  problem  reduces  to 
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(E-10) 


r 


d 

2 f 

2g 

0 

0 

2f 

2k 

2 f 

2h 

0 

2g 

2 f 

2d 

2 f 

2 o 

0 

2h 

2 f 

2k 

2 f 

0 

0 

2g 

2f 

d 

2B? 

2153 

215, 


The  matrix  on  the  left  side  is  decomposed  via  a Choleski  decomposi- 
tion, and  the  resulting  equation  is  solved  using  the  procedure  outlined 
in  Eqns  (4-6)  and  (4-7).  Also,  with  the  definition  of  T as  indicated  in 
E < i (E-3),  the  values  of  T calculated  from  Eqn  (E-lO)are  actually  the 
values  of  the  potential,  cj> . 


Decoup 1 ing 

An  attempt  was  made  to  reduce  the  large  matrix  problem  to  two  smaller 
problems  by  "decoupling".  If  Eqn  (4-llb)  is  used  to  approximate  the  angular 
derivatives  at  each  mesh  point,  every  other  radial  line  is  linked  together. 
All  the  mesh  points  along  the  even-numbered  radials  are  linked  together, 
and  all  the  mesh  points  along  the  odd-numbered  radials  are  linked  together. 
Thus,  for  the  sample  mesh  in  Fig  E-l,  the  40  x 40  matrix  could  be  reduced 
to  two  20  x 20  matrices.  (These  could  be  further  reduced  through  the  use 
of  symmetry. ) 

This  technique  did  not  work  because,  as  explained  in  Chapter  IV,  use 
of  Eqn  (4-llb)  is  an  improper  way  of  approximating  the  angular  derivatives. 
Bo  left  and  right  derivatives  must  be  approximated  at  each  mesh  point, 
and  this  destroys  the  "decoupling"  nature  of  the  matrix. 
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Appendix  F 


Importance  of  Na tu rnl  Boundary 


Conditions  when  Minimizing,  a Functional 


The  first  attempt  at  the  variational  approach  was  to  minimize  the 


funct ional 


s = j f(v&yciA 

"a 


(F-l) 


and  impose  the  Neumann  boundary  conditions  through  the  use  of  a Lagrange 
Multiplier,  A. 

In  matrix  form,  the  attempt  was  to  minimize 


I T 


6 A 


( F-2) 


subject  to  the  constraint 


a = o 


( F-3) 


where  B and  a were  determined  from  the  quadrature  formulae  and  the  boundary 
conditions . 

Imposing  Eqn  (F-3)  on  Eqn  (F-2)  through  the  use  of  Lagrange  Multipliers 
led  to  a new  functional,  S*: 


s*  - 2 A ± + a (8  £ - a) 


(F-4) 


(Note:  A is  used  to  be  consistent  with  the  definition  of  matrix  multiplica- 

tion. ) 
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Then 


S*  was  varied  with  respect  to  <J>  and  X. 


- 4^  + J*ii  - ° 


a 6$  - a - o 

i >1  — 


(F-3) 


(F-6) 


Matrix  Eqns  (F-5)  and  (F-6)  were  then  solved  to  determine  $ . Figure  F-l 
illustrates  the  results  of  the  calculation. 

The  error  in  the  above  procedure  is  that  the  functional  in  Eqn  (F-l) 
has  inherent  in  it  homogeneous  Neumann  boundary  conditions.  This  can  be 
seen  by  setting  g0  equal  to  zero  in  Eqn  (D-3).  Thus  the  functional  was 
trying  to  minimize  an  equation  with  inherent  (natural")  homogeneous  Neumann 
boundary  conditions  on  both  boundaries.  But,  through  the  use  of  Lagrange 
Mu  ’ipliers,  it  was  also  being  told  to  satisfy  inhomogeneous  Neumann 
conditions  on  the  outer  boundary.  These  two  conflicting  conditions  pro- 
duced the  unsatisfactory  curve  in  Fig  F-l. 

An  investigation  was  also  made  to  determine  the  effect  of  imposing 
Neumann  boundary  with  Lagrange  Multipliers  on  a functional  in  which  the 
Neumann  conditions  were  already  imposed  through  the  natural  boundary  condi- 
tions. The  following  one  dimensional  problem  was  used  for  this  purpose: 


R'ta)  = o 

R'W  = 1 


The  corresponding  functional  is 


( F— 7 ) 
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7 DO  ~ 


In  matrix  notation,  the 

/A  R.  "?*  C A JB  (F-9a) 

CR  - a =:  O (F-9b) 

where  A = coefficient  matrix  from  the  integral  term 

C = coefficient  matrix  generated  from  the  boundary  conditions 
13  = vector  generated  by  the  second  torin  of  the  functional 
a.  = boundary  condition  vector 

As  illustrated  in  Fig  F-2,  the  additional  enforcement  of  the  boundary 
conditions  through  the  use  of  Lagrange  Multipliers  magnifies  the  scallops. 
For  this  reason  the  two-dimensional  problem  was  solved  by  enforcing  the 
Neumann  conditions  through  the  use  of  natural  boundary  conditions. 


i( 

r 


T 


(rT  -i-  £ 


2 RM  0) 


rd' • + 


( F-8  ) 


problem  is 


2 3 4 

Radial  Distance 


F.ig  F-l.  Potentials  alcr.g  the  Line  0 = 0 as  Solved  by 

Imposing  the  Boundary  Conditions  with  Lagrange 
Multipliers  and  Ignoring  the  Natural  Boundary 
Conditions. 
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1 


2 

Radial  Distance 


3 


Fig  F-2.  Solution  of  the  ,1-D  Problem  Using  Lagrange 

Multipliers  and  Natural  Boundary  Conditions. 
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Append  Lx  G 

Variational  Computer  Program 

The  purpose  of  this  appendix  is  to  explain  the  computer  program  used 
to  solve  the  matrix  equation  generated  by  the  approach  using  variational 
calculus.  The  main  program  and  each  subroutine  will  be  explained  in  a 
separate  sub-section.  The  format  includes  an  explanation  of  all  variables, 
references  to  all  equations,  and,  when  necessary,  a card-by-card  explana- 
tion of  the  program. 

Main  Program 

The  purpose  of  the  main  program  is  to  call  various  subroutines  which 
first  set-up  and  then  solve  the  matrix  equation. 

Q = Coefficient  matrix  (stored  as  a linear  array). 

Y = Boundary  condition  vector. 

NR  = Number  of  radial  mesh  points;  must  be  an  odd  integer. 

NA  = Number  of  angular  mesh  points;  must  be  an  even  integer. 

RI  = Inner  radius  of  the  cylinder. 

RO  = Outer  radius  of  the  cylinder. 

Cl  = Incoming  current  density. 

H = Radial  mesh  spacing. 

W = Angular  mesh  spacing. 

NQDIAG  = Length  of  the  main  diagonal  of  the  coefficient  matrix.  It  is  the 
minimum  allowable  dimension  for  the  Y array. 

NQL  = Minimum  allowable  dimension  for  the  Q array. 

SIGMA  = Conductivity  of  the  cylinder. 
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Cards  4,  5:  LLI.  is  any  non-zero  integer.  Its  purpose  is  to  allow  for  the 
input  of  more  than  one  set  of  data.  A card  containing  a value  for  LLL 
must  be  plciced  before  every  card  containing  input  data.  The  last  card 
in  the  data  deck  must  be  the  integer  zero.  This  will  stop  the  program. 

Cards  13,  14:  The  values  of  the  potential  are  printed  out. 

Subroutine  Data 

This  subroutine  reads  in  the  data  and  calculates  all  the  parameters 

required  for  the  solution  of  the  matrix  problem. 

Card  2:  Input  data  is  read  in.  Free  format  is  used. 

Card  5:  The  length  of  the  Y vector  is  calculated.  This  number  is  the  minimum 
allowable  dimension  for  array  Y. 

Card  6:  Q is  a symmetric  matrix.  The  lower  triangular  of  Q is  stored  in  a 
linear  array,  the  length  of  which  is  determined  by  the  given  formula.  It 
is  the  minimum  allowable  dimension  for  array  Q. 

Subroutine  QSETUP 

This  subroutine  sets  up  the  coefficient  matrix,  Q. 

Cards  3,  4:  The  Q array  is  zeroed. 

Card  5:  A DO  LOOP  is  initiated  for  the  calculation  of  the  coefficients  of 
the  mesh  points . 

Card  6:  A subroutine  is  called  which  calculates  the  appropriate  Simpson's 
Rule  coefficient  for  each  mesh  point. 

Card  7:  The  radial  distance  to  each  mesh  point  is  calculated. 

Cards  8,  9:  The  constant  coefficients  from  Eqns  (4-9)  and  (4-11)  are  calculated. 

Cards  10  - 20:  The  appropriate  subroutine  for  is  called  for  the  mesh  point 
under  consideration,  depending  on  whether  it  lies  along  the  inner  radius, 


the  outer  radius,  or  somewhere  in  between. 


Card  22:  Since  the  matrix  is  singular,  the  lower  left  element  will  be 
zero  after  the  Choleski  decomposition  is  applied.  This  card  sets  that 
element  to  zero. 

Sk ibroutinc  SC 

This  subroutine  calculates  the  Simpson's  Rule  coefficient  for  eacli  mesh 
po int . 

Cards  2-10:  Modular  arithmetic  is  used  to  determine  the  radial  position 
of  each  mesh  point,  and  assign  a value  according  to  the  following  pattern: 

1,4,2,4,.. .4,2,4, 1. 

Card  11:  Since  a two-dimensional  mesh  is  being  used,  Simpson's  Rule  has  to 
be  applied  in  the  angular  direction  also.  This  card  assigns  a value  of 
2 to  mesh  points  along  odd  radials,  and  a value  of  4 to  those  lying  along 
even  radials. 

Card  12:  The  angular  and  radial  contributions  to  the  Simpson's  Rule  coef- 
ficient are  multiplied  together. 

Subroutine.  QINNER,  QOUTER,  QMIDDLE , and  QJUNC. 

These  subroutines  calculate  the  contributions  of  each  mesh  point  to  the 
\ Q matrix.  Eqns  (4-10)  and  (4-11)  are  used  to  generate  the  equations  in  the 

subroutines.  The  constant  term  was  not  included  in  these  subroutines  because 
it  was  accounted  for  in  the  calculation  of  Cl  and  C2  in  QSETUP.  Since  linear 
symmetric  storage  is  being  used,  it  is  necessary  to  convert  from  a two-dimen- 
sional array  to  a linear  array  using  Q (l,J)  = Q(K)  = Q(I  ( I— 1 ) / 2 + J). 

Subroutine  LSETUP 

This  subroutine  performs  a Choleski  decomposition  of  the  Q matrix.  The 
diagonal  elements  are  calculated  using 
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(0-1) 


The  remaining  elements  are  calculated  using 


4*  iN -£'«'») 


(0-2) 


Cards  3,  5:  A DO  LOOP  is  initialized  which  spans  all  but  the  final  element 
of  the  lower  triangular  matrix.  This  element  has  been  set  to  zero  in  t he 
previous  subroutine. 

Card  6:  111  is  used  to  span  the  elements  in  the  first  column  of  the  matrix. 

Card  7:  JB  is  used  to  span  the  row  elements. 

Card  8:  ID  is  used  to  span  Che  main  diagonal  elements. 

Card  9:  11  is  a parameter  used  by  VPROD  to  indicate  the  number  of  pairs  of 

elements  being  multiplied  together.  It  is  always  one  less  than  the  column 
index  of  the  element  being  calculated. 

Card  10:  Eqn  (G-l)  is  used  to  calculate  the  main  diagonal  elements. 

Cards  LI,  12:  A DO  LOOP  is  initialized  which  spans  all  elements  beneath  the 
main  diagonal  element  just  calculated. 

Cards  13  - 15:  Eqn  (G-2)  is  used  to  calculate  the  elements  below  the  main 
diagonal . 


S u brou tine  YSET 

This  subroutine  calculates  the  inhomogeneous  boundary  condition  vector. 
It  is  the  right  side  of  Eqn  (4-5).  (Its  infinite  sum  representation  is 
shown  in  Eqn  (4-13).)  The  only  non-zero  values  of  this  vector  will  be  those 
corresponding  to  a mesh  point  lying  on  the  outer  radius. 


66 


1 

Cards  3,  4:  The  Y array  is  reset  to  zero. 

Card  5:  A DO  LOOP  is  initialized  which  spans  all  mesh  points. 

Card  6:  II  the  mesh  point  under  consideration  does  not  lie  along  the  outer 
boundary,  the  DO  LOOP  continues  without  assigning  it  a value. 

Cards  7,  8:  The  appropriate  Simpson's  Rule  coefficient  is  determined 
depending  on  whether  the  point  lies  along  an  odd  or  even  radial. 

Card  9,  10:  Eqn  (4-18)  is  used  to  calculate  the  boundary  value.  In  this 
case , f ( 0 ) - cos  0 . 

Card  12:  The  last  value  of  the  vector  is  arbitrarily  set  to  zero,  as 
explained  in  the  last  paragraph  on  page  14. 

Subroutine  SOLVE 

This  subroutine  solves  the  matrix  equation  L L = % by  first  calculating 

x in  Lx  = y using  forward  substitution,  and  then  calculating  ij>  in  L^73  x using 

backward  substitution. 

Card  5:  A DO  LOOP  is  initialized  which  spans  all  but  the  lower  left  element 
of  the  matrix. 

Card  6:  11  is  used  with  VPROD  to  indicate  the  number  of  pairs  of  elements 

being  multiplied  together. 

Card  7:  IB  is  an  index  which  enables  VPROD  to  calculate  the  row  products  of  L. 

Card  8:  ID  spans  each  element  of  the  main  diagonal. 

Card  9:  The  values  of  x are  calculated  and  stored  in  the  Y array. 

Card  10:  The  last  element  of  x is  set  to  zero. 

Card  11:  This  card  initiates  the  back  substitution  scheme  by  calculating  the 

second  last  element  of  £ and  storing  it  in  the  Y array. 

Cards  12,  13:  Since  the  last  element  of  the  £ vector  is  zero,  and  the  second 
last  element  has  just  been  calculated,  the  DO  LOOP  must  be  indexed  to  span 
NQDIAC-2  elements. 
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Card  14:  I is  used  to  index  the  calculation  of  the  Lth  element,  of  <p . 

Card  16:  L spans  each  element  of  the  main  diagonal. 

Cards  17  - 19:  This  loop  calculates  the  column  products  of  L with  the  cor- 
responding terms  of  the  $ vector. 

Card  20:  The  terms  of  the  ^ vector  are  computed  and  stored  in  the  y array. 
Cards  21  - 26:  The  vector  is  normalized  as  explained  in  the  last  paragraph 
on  page  14. 
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Appendix  H 


■ 'er  i : 1 tj  on  _o_f_  Che  Green's  Fu notion 
and  Iteration  Formulas 


This  appendix  contains  three  derivations.  First,  the  form  of  the 

V -> 

Green's  function,  G (x-x1)  is  derived.  Secondly,  detailed  iterative 
fori.ulas  art?  derived  for  i|/(x)  and  <ji(x)  as  used  in  Eqns  (5-8)  and  (5-9) 
respect ively . 

D e _ri v . it:  ion  o f t h e Gr eon's  F 1 1 n c t ion 

As  stated  in  Chapter  V,  if  G (x-x1)  is  a solution  to  V-(J>  = 0,  then 

V' Gl  (X - k ■ ' ) a.  S (‘x  - “£')  < H - D 


Since  the  problem  is  being  solved  in  cylindrical  coordinates,  the  delta 
function  takes  the  form 


cf(r) 
-Tl*  T 


I -v  -y  . 

where  r = | x-x  | 

This  equality  can  easily  be  seen  from  the  following: 


( H — 2 ) 


(11-3) 

j 

where  f(x)  is  a test  function,  (Ref:  11,  2.9). 

In  terras  of  an  integral  over  r and  6,  the  left  side  of  F,qn  (H-3) 
becomes 
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l?r  ro 

J J-f(r)  r dr  do 

o o 

? r 

~ 7 (O) 

an 

zj fCO  cf  <*)  dr 

O 

z[i  *(0\ 

= -H°) 

The  equality  holds,  therefore  Eqn  (ii-2)  is  valid. 

Substituting  Eqn  (H-2)  into  Eqn  ( I! — 1 ) , and  writing  in  cv 
coordinates  yields 


.L  (r  - \ . cf 

7 sv  \ a r / *r  r 2 ai-1  ~ i? 


t a Ci 


cf£0 

r 


Assuming  the  Green's  function  does  not  depend  on  angle,  the  part 
tive  with  respect  to  0 vanishes.  Thus  Eqn  (li-7)  becomes 


Integration  of  Eqn  (H-8)  yields 


5.L  f r £!& ) ~ d'  r ) 

d r \ ci  V / it* 


y»  .d  '~~i  — 
c;  r " ir 


Equation  ( 11— 9 ) is  easily  solved  for  G: 


where 


-v-fc  m 

Gs  U-x  ) = 


vV.x.x’)  = o 


In  i* 


1 

( H—  -i- ) 

(11-5) 

(11-6) 

1 i nd  r ica ) 

(H-7 ) 

i.al  deriva- 

(H-8) 

(H-9) 

(H-10) 
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r - |«-  x'| 

KuC  (3,  92-102)  contains  a detailed  discussion  of  logarithmic  potentials, 
including  derivations  which  amplify  on  the  factor  of  x/i\. 

If  H(.x,.x')  is  set  to  zero,  the  gradient  of  G is 


VO,  (x-  *')  = V 


* (.±\ 

M'r  / 


( H — 11) 


(H-12) 


C H— 13) 


Since  the  Green's  function  is  symmetric 


V '<a  ( S'-  -ILzJL. 


01-14) 


Equations  (H  10)  and  (H-14)  are  used  in  conjunction  with  Eqns  (5-8) 
and  (5-9)  to  solve  for  d> . 


Iterative  Formula for  iji(x) 

The  formula  for  ij>  is  given  hi  Eqn  (5-8)  as 


V C?  ) a -Jg  (?-?)  V <j> (*>  d a' 


(11-15) 
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Substituting  Eqn  C H— 10)  and  per  topping  the  dot  product  operation 
y ie Ids 


V(x)  ^ 


f In  lx '-3(| 


'J 

S 


ir 


' Q 

— a s 
2X 


(H-16) 


Since  there  is  only  a contribution  from  the  outer  boundary,  Equ  (H-16) 


can  be  simplified  to 


( H — 1 7 ) 


Or,  as  a summation  over  N mesh  points  on  the  outer  boundary 


N 


^(x)  In  i 


n = < 


(H-18 ) 


When  calculating  T on  the  inner  boundary  due  to  the  source  points  on  the 
outer  boundary,  the  above  formula  presents  no  problem.  However,  when 
calculating  'F  on  the  outer  boundary,  the  term  |x’ -xj  approaches  zero  as 

| — V 

x approaches  x.  Fortunately  this  "self-contribution"  term  can  be  calcu- 
lated analytically. 

Figure  (H-l)  illustrates  the  schematics  for  the  calculation  of  the 
self-contribution"  term.  The  point  iji  cannot  be  evaluated  numerically. 
Points  £ and  <j>_  are  neighboring  mesh  points. 

If  the  points  are  close  together,  0'  is  very  small,  and  3<|>/3r  is 
approximately  constant  over  the  interval.  So  for  a small  interval, 
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8! 


r ’ 


ihis  express  i.on  must  b added  to  the  numerical,  non-n  iuigular  summation  of 
Eqn  (K-L8). 

Iter  ) t tv_*  i J*\  > r t n ; !.  a for  <fi  ( x ) 

Tiie  i o rmu  1. a for  $(x)  is  given  in  Eqn  (‘5—9).  If  the  surface  integral 
ovol  S is  factored  into  Integra’s  over  the  inner  and  outer  boundaries, 
the  equation  can  be  written  as 


m ' » 4 \ 


-r-“  1 <• 

* * i V / 

) rK  s\ 


4* 


(H-22) 


r_\  <fc5>  t . 

J ft'-W 

b 

where  Eqn  (H-14)  has  been  subst L tutod  for  VG. 

Bor  this  derivation,  the  values  of  <•)  will  be  evaluated  on  the  inner 
boundary.  figure  H-2  illustrates  the  method  of  evaluating  the  integrals 
in  Eqn  (H-22).  Consider  first  the  contribution  from  all  points  on  the  inner 
boundary.  Once  tne  dot  product  operation  has  been  performed,  the  first  inte- 
gral of  Eqn  (H-22)  can  be  written  as 


r A / V ')  £ , /Tf*  „ -y  ^ 

/ a •) 

J V j ;;  - ;;  | 


(H-23) 


ct 


As  can  be  seen  in  Fig  H-2,  the  angle  between  x^'—x  and  n is  — ot.  Also  since 
the  triangle  O-x-'-x  is  isosceles,  2a+0=-r.  Thus 


A -a  / — !> . ^|  | 


XQ-  x j cosf-jf-) 


32. 


(H-24) 


(H-2 5) 


(11-26 ) 


1 ''a.  ■'  J ?.  ; 


Further,  utilizing  the  Law  of  Cosines 


\K^x\  =.  z o?  - z c»*  o 


( 11—2  7 ) 


After  factoring  23^  and  applying  a half-angle  identity  for  the  sine, 


i >?«  * ■ * i ~ ^ a? 


(H-28) 


When  Eqns  (11-26)  and  (h-28)  arc  substituted  into  Eqn  (H-23)  the  result  i? 


r.*&)  fjs 

5?  rr 


J 

a. 


f O'JMds 


(H  - 29) 


However,  as  a consequence  of  the  boundary  condition 


f-tds  ~ 0 


J dr 


(H-30) 


Eqn  (11-29)  is  also  zero.  Hence  it  is  not  necessary  to  evaluate  this  terra 
when  calculating  <j>. 

Next,  the  contribution  from  all  points  on  the  outer  boundary  is 
determined.  The  second  integral  of  Eqn  (11—22)  can  be  written 


r <h?‘)  ! tj-  *1  c-°s  & c\  3 


(11-31 ) 


• | -r 

where  8 is  the  angle  between  -x  and  the  outward  normal  at  the  source 
point.  The  above  integral  car.  be  written  as  the  following  summation, 


b ?.‘'j  V~*  C.PS 
~ i~.  jx'-xf 

n = i * t>  ■ v\ 


(H-32) 


a 
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where  8^0 nd  x^'-x  can  be  detoi mined  Lcci  the  Law  of  Cosines. 

liie  iteration  formula  lor  ij>(;<)  on  the  innc-r  boundary  can  new  bt 
found  by  substituting  Ecjn  (11-32)  into  Eqn  (11-22). 


PJ 


■.o-'  p <t>  (:zy  a 


,S  /v'\  - <>  f y \ • b f. 

fcn 


(11-33) 


A ill  lilar  formula  for  ealculat  ing  ij)(x)  cii  the  outer  boundary  can  be 
found  by  following  a similar  procedure.  The  formula  is 

M 


<k  (?)  = 

o 

VJ  it) 

‘fc  ' ' 

1 'Tr  / 1 c/ 

•*  C i i ^ ' 

V 1 

3 >4  | 

ns/  i -y 

1 r\ 

where  y is  the  angle 

between 

-y 

>'-b  arid  the  outward  unit  normal  at  the  source 

points . 


A computer  program 
iterate  values  of  ip  can 


using  Eqns  (11-18),  (H-21), 
be  found  in  Appendix  J. 


(11-33),  and  (H-34)  to 
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AD-A064  754  AIR  FORCE  INST  OF  TECH  WRIGHT-PATTERSON  AFB  OHIO  SCH— ETC  F/G  20/3 
SPACECRAFT  CHARGING. (U) 

DEC  78  R A PASSOW 


UNCLASSIFIED  AFIT/GNE/PH/78-19  NL 


As  6 approach":!  ze  ro,  i he  integral  on  Che  right  approaches  .>(x) , :o 

upon  substitution  of  Eqns  (1-7)  md  ( I — 1 0 ) into  Eqn  ( I— 6 ) 


r. 


j 


In  Iy-'*!'?3#  <M  - 


[ln|y-H|  It  - in  Iv-*l 


cJS  ■'*.  £ tr i*)  (i-n) 


Or, 


${*)  =■  A;  fin  ly-xj  — 


J 

■A 


y% 

\;j  In  i y- X i ||  - 4>|^  In 


(1-12) 


Since  = 0, 


K*)  = 2Wr*<y)|i^lv-'4-!nlv-.v!||]lr.  (*-»» 


Consider  a point,  n,  on  the  boundary,  S.  As  x approaches  q 


lim 

x->*v 


* A r v'“  Sn  i n I y ” ! d S = 


Tr<!Hs0  +J^(V)  ~ In  1 y-n't  ds 


(1-14) 


Thus,  Eqn  (1-13)  becomes 


r 

4>(n)  = 4 4><wi  i-  <*><*&  •"  ,y'nl  ds  “ 


2-tfj 
.*? 


l 

2 IT 


:jln!  y-hiH^S 


(1-15) 
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Note:  Eqn  (L-14)  is  derived  in  Appev  !ix  A of  Reference  3. 
Fi nal ly , 

4>(n)  = ln|y-nl  ds  - 

S 

^ |y-nlH  ds 

s 

Eqn  (1-16)  is  valid  at  the.  boundary,  and  if 


(1-16) 


x s n 

smh  f 

X a y 


(1-17) 


(1-18) 


G s in  iS'-xl 


(1-19) 


it  is  seen  to  be  the  same  as  Eqn  (5-7). 


1 


Appendix  J 

Green 1 s P u ' i ction  Comp u ter  Progra m 

This  program  is  divided  into  seven  main  parts  by  groups  of  three 
comment  cards.  Basically,  the  program  follows  the  iteration  scheme  out- 
lined in  Chapter  V. 


Part 

1: 

Input  data 

i.s 

read 

in  and  all  required  constants  are  calculated. 

Part 

2 ; 

The 

value 

of  ‘p  on  th 

e inner  and  boundaries 

is  ca 

lc  ilated  as  given 

in  Eqns  (H 

-18) 

and 

(H-21)  . 

Part 

3 : 

Eqn 

(11-33) 

is 

used 

uo  calculate  $ at  each 

i nni  r 

bounlary  mesh  point. 

Part 

4: 

Eqn 

(P-34) 

is 

used 

to  calculate  <j>  at  each 

out«'r 

boundary  mesh  point. 

Part 

5: 

All 

values 

of 

<J>  are 

tested  for  convergence. 

Part 

6: 

Eqn 

(5-14) 

is 

used 

to  calculate  <j>  in  Region  II. 

Part 

7: 

The 

boundary  conditions  are  calculated  as 

gi  ven 

in  Eqn  (2-21). 

( 
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FPOG!-  AM  GPEEHC  INP'.J  r , 0(JTPUT> 

DIMENSION  SA('!T]),Gl(?ro),PA(20J),P3(200), 
1 PAS'AV'  (?H0)  ,3IBSA  VP(2G0) 

READ1  , A, A ,01  , SIGMA , M A 

SS  =A  •*  A (■■P  0 

ST  = 2 .*  A*3 

S3S  = 2 ,-'8*3 

SAS=?  A 

n 3 = A » n 

AS=A* A 

PI  = 3 , It,  lh>«2oS.3S) 

DA=2.*FI/MA 


CA  LC  PLATE  PSI  ON  I 'HER  A NO  OUTER  BOUNDARIES. 
HO  1 1=1, N A 
S ° ( I ) = 0 . 

SA (I) =0. 

CO  3 J = l » N A 

<S  , - 2 ,»M30(  Jvl,  ? ) ) / 3 . 

THETA  = (T-J)’DA 

SA  (I)  = SA  ( T)  fSIMp  ' P'.n  C i,OA  ,01  ,SIO-<A)  » 

1 ALOG (SOcT ( S ? - S r " C G S ( THETA) > I 

IE(I.EO.J)GO  TO  3 
I E { M ( 0(1  f.J,2)  . EO  .Cl)SrMP=A  ,/3  * 

IF(M(D(I+J,c)  .NE.(l)SIMP=2./3. 

TE(I  A RS(  -00(1,  NA)  - ‘in  ( J,NA)  ) . EO.  1.  DR. 

1 I A 0 S (VOO( I, vA ) -MOP ( J,MA ) ) 

2 . E0.NA-1)SIM3  = 1./7. 

SP(I)  ~ S ’■>  ( I)  *ST  M°‘  OHO  ( J * D A , C I , SIGMA)  * 

1 ALOG (SORT (SES -SOS'  COS( THETA) ) ) 

3 CONTINUE 

S3  (I ) =-8*'  0A/Pta53<  I) 

SA (I) =-E>DA/PI*SA ( I) 

S3 (I ) = S'J ( I) ~2.  f 1NO (I , DA , Cl, SIGMA) /PI  * 
l P* DA* (ALOG (P* OA) -1) 

1 CONTINUE 

DO  L 1=1, NA 
PA  SAVE ( I ) =SA ( I ) 

°3SA  V T ( I ) =SQ ( I ) 

PA (I ) =SA ( I) 

4 P 3 ( I > = S 0 ( I ) 
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CALCULATE  °HI  ON  INNER  BOUNDARY. 
rQ  5 1=1, N A 
T E U P n = 0 . 

6 J = i » M A 

31 (<-.-?.»  >K)0  ( j *■  T , ?))/3. 

THETA = (I-J) “DA 
XS=S1 -ST* COS (THETA ) 

COSA = (VS*DS-AS)/(? .‘l* PORT  (XS) ) 

T EMP  '■  = rENPDtSI  MO*  SOS  A "PR  ( J)  /SORT  ( <S) 
CONTIf HE 

PA  C)  =SA  (I)  ♦•’3* QA/P I1* T F PPQ 
CONTINUE 


CALCULATE  -HI  ON  OMTr-R  -BOUNDARY. 

PC  11  I = 1 1 N A 
TEUPA=0. 

TO  11  J = 1 , NA 

SIHO  = (<- J+I,2)  >/3. 

THET t = ( C - J ) ‘ 0 A 
XS=Sf-ST ' COS ( THETA) 

COSD= (YSvAS-OS) / (0 , * A * c 0 R T ( X S ) ) 
C0SGrCGS(°I-ASCSC0S0)  ) 
TEMDArTERPAtSIH0*  SOSO*  PA  (J) /SORTUS) 
CONTINUE 

P'’  ( I ) = S 3 ( I ) +A*  T FNPA 

CONTINUE 
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TEST  ;rCr<  CONVERGENCE. 

21  1=1  , NA 

TESTA  = APS  (PA  (I ) --•'ASIVE(I)  J 
TEST”  = CPS  ( " H (I  ) - - 3 S \ V E ( I ) ) 

IF(TC'SI"\.GT. .30  9 15  GO  TO  21 
2 C IE (TFST9  .GT.  .33  31)  GO  TO  21 
GO  ro  22 

21  DO  2?  1=1 , N A 

°CSAVl ( I ) =PA (I) 

2?  DRSAVE(I) =DD(I) 

PRINT* , "ONE  ITERATION'* 

GO  TO  2 

2?  DO  2 <■  1=1  , NA 

?s  PRINT*  ,°A  (I)  ,PP(I>  ,**  M A ='* , NA 

C 
C 

C CALCULATE  PHI  IN  REGION  II. 

READ*  ,PP 
no  31  J=  1 , N A 
THETAP=(  J-l) '‘360,/NA 
TEHP=  0 

DO  3 C J = 1,NA 

TP  ETC  = (T  - l)  * c.  /NA 

TTP=  (THEIAP-TH  ET  !)»?.*  PI/ 3S  0 . 

9RPS  = RF  + RD«-AM-2.  “R351  A-'  COS (T TP) 

C 0 £ G = (RP3St£S-RP'  RR) / < ?.*A* SORT (RRPS)  > 
Tr  TENPsTE.M3*C0SG-'Of  ( I)  /S^RT  (RRRS) 

TEU3=T EMP*  A * OA/°I 
T 1 PRINT*  ,RP,THETAP,rENP 
END 


rjNOIIGN  °NO(J,Of ,CI,STGMA) 
CT=CCS ( ( J-l) »DA) 

IF  (C  .LT  . 0 .>  CT=3  . 

r*ND=F  I /GIG'^'A  ' (CT-.  31  8309835164) 

P E TUT N 

ENO 
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Richard  Alan  Passow  was  born  in  Milwaukee,  Wisconsin  on  19  June  1950. 
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The  spacecraft  is  modeled  as  a hollow,  infinite  cylinder.  The  fields 
are  calculated  for  the  case  in  which  the  incoming  current  is  incident  per- 
pendicularly to  the  longitudinal  axis  of  the  cylinder.  Basic  electro-static 
theory  reveals  that  the  governing  equation  for  the  potential  is  Laplace's 
Equation,  subject  to  Neumann  boundary  conditions.  This  equation  is  first 
solved  by  separation  of  variables.  The  Ii-Field  predicted  for  representive 
values  of  incoming  current  density  are  on  the  order  of  IQ-7  volts/meter. 

Several  pitfalls  encountered  with  the  variational  approach  are  explained. 
These  include  the  importance  of  natural  boundary  conditions,  ensuring  continu- 
ity of  the  first  derivative  across  all  mesh  points,  and  difficulties  encoun- 
tered in  trying  to  reduce  the  size  of  the  matrix  through  "decoupling".  The 
results  of  this  section  show  that  the  variational  method  tends  to  approximate 
the  analytic  solution  as  the  mesh  becomes  finer. 

9 

The  Green's  function  approximations  of  the  potential  distribution  were 
consistent  with  analytic  results  to  at  least  two  significant  figures. 
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